//Luokasta cNum peritty cRK4 - Runge-Kutta - menetelmän toteutus.
//
//Joonas Koskinen
//Pro Gradu Työ
//Jyväskylän Yliopisto

// 	0       |
//	1/2     |        1/2
//	1/2     |        0 	1/2
//	1       |        0 	0 	1
//      ------------------------------------
//              |        1/6 	1/3 	1/3 	1/6

#include "cRK4.h"

RK4::RK4() {
	//double ivy, ivdy;
	setName("RK4");
	//tx and fx are cleared and set from cInfo when active potential is set.
	//If, for some reason, the active potential is not set and method is initiated we 
	//set tx and fx to 0 and 10 respectivly.
	//ivy=0.1;
	//ivdy=0;
	//setStep(0.01);
        /*cout << "RKN set step to " << getStep() << endl;
	set_tx(0);
	set_fx(10);
	set_ty(ivy);
	set_fy(ivy);
	set_tdy(ivdy);
	set_fdy(ivdy);*/
}

RK4::RK4(double ivy, double ivdy) {
	setName("RK4");
	//tx and fx are cleared and set from cInfo when active potential is set.
	//If, for some reason, the active potential is not set and method is initiated we 
	//set tx and fx to 0 and 10 respectivly.
	setStep(0.01);
	set_tx(0);
	set_fx(10);
	set_ty(ivy);
	set_fy(ivy);
	set_tdy(ivdy);
	set_fdy(ivdy);
}

RK4::~RK4(){}

void RK4::methodStep(double& step, vector<double>& xx, vector<double>& yy, vector<double>& dy, cDY* DY, cPot* V) {
	double k1,k2,k3,k4,l1,l2,l3,l4;
	k1 = dy.back();
	l1 = DY->DE(xx.back(), V->potential(xx.back()),  yy.back(), dy.back());
	k2 = dy.back() + step * 0.5 * l1;
	l2 = DY->DE(xx.back() + 0.5 * step, V->potential(xx.back()+0.5*step), yy.back()+0.5*step*k1, dy.back()+0.5*step*l1);
	k3 = dy.back() + 0.5 * step * l2;
	l3 = DY->DE(xx.back() + 0.5 * step, V->potential(xx.back()+0.5*step), yy.back()+0.5*step*k2, dy.back()+0.5*step*l2);
	k4 = dy.back() + step*l3;
	l4 = DY->DE(xx.back() + 0.5 * step, V->potential(xx.back()+step), yy.back()+step*k3, dy.back()+step*l3);

	yy.push_back(yy.back() + step*(k1 + 2.0 * k2 + 2.0 * k3 + k4)/6.0);
	dy.push_back(dy.back() + step*(l1 + 2.0 * l2 + 2.0 * l3 + l4)/6.0);
	xx.push_back(xx.back() + step);
}
